Numerical Modeling of Charged Black Holes 
with Massive Dilaton * 

T.L. Boyadjievj P.P. Fiziev* 



Abstract 

In this paper the static, spherically symmetric and electrically charged black hole solutions in 
Einstein-Born-Infeld gravity with massive dilaton are investigated numerically. 

The Continuous Analog of Newton Method (CANM) is used to solve the corresponding non- 
linear multipoint boundary value problems (BVPs). The linearized BVPs are solved numerically 
by means of collocation scheme of fourth order. 

A special class of solutions are the extremal ones. We show that the extremal horizons within 
the framework of the model satisfy some nonlinear system of algebraic equations. Depending on 
the charge q and dilaton mass 7, the black holes can have no more than three horizons. This allows 
us to construct some Hermite polynomial of third order. Its real roots describe the number, the 
type and other characteristics of the horizons. 

Keywords: black hole, scalar-tensor theory of gravity, massive dilaton field, multipoint boundary 
value problems, continuous analog of Newton method, method of collocation. 
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1 Introduction 

Being a product of strong nonlincarity of the modern theory of gravitation, the black holes (BH) 
represent a significant challenge since their theoretical discovery within the framework of a general 
theory of relativity (GTR) p^. Till now BH are not a subject of a direct experimental analysis and 
at present there is no indisputable example of observed BH with their critical property - the existence 
of an event horizon 0. Moreover, there exist serious physical arguments supporting Einstein-Dirac's 
assertion 3, that the horizons are not a physically admissible notion. In GTR there also exist a 
large number of solutions of the corresponding physical problems without horizons at all. A more 
detailed discussion of these problems can be found in 4 . In this complicated situation one needs a 
trustworthy methods of both theoretical and experimental study of BH to make a serious and well 
founded conclusions about the real physical meaning of the solutions with horizons and their ability 
to describe the physical reality. 

The modern theories of fundamental interactions, such as the theories of space-time with torsion, 
dilaton gravity, Born- Infold electrodynamics, supergravity, superstrings, superbrane and M-thcory have 
multiply enriched the set of different BH solutions. These models have not yet reach their final form 
and do not have a status of phenomenologically justified physical theory, nevertheless, they present 
some theoretical interest at least as a mathematical examples of nonlinear theories. 

From physical point of view BH can be divided into two basic classes: 

1. Macroscopic BH, with the mass varying from Chandrasekhar's mass ~ 3M Q to mass ~ 10 6 — 
10 11 Mq, i.e., about mass of a Galaxy; 
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2. Microscopical BH, whose mass has the order of that of elementary particles of the Standard 
model. 

If today there are some preliminary indications that the macroscopic BH may be located in active 
cores of galaxies, in the (not active) core of our Galaxy, in quasars, and also as a dark components in the 
systems of X-ray star pairs [5], then for the time being we have not at disposal any phenomenological 
indications about an existence of microscopical BH. Until recently the basic reason for that was the 
circumstance, that the microscopical BH represent a non-linear formation in 4D space-time, whose size 
is of order of the Planck length ~ 10~ 33 cm. This is far outside of limits of experimental accessibility 
even in distant future. 

During the last years the situation in that plan has varied a little: miscellaneous theoretical models 
in the frameworks of superstrings and superbranes augmented essentially the anticipated size of possible 
microscopical BH at the expense of usage of (at present not observed) higher dimensions of space- 
time. Therefore now one considers the opportunities for discovery of microscopical BH on the future 
accelerators like LHC, VLHC, NLC at energies about several TeV (see, for example, E] and the 
references therein), and the generation of such kind BH in the cosmic rays also (see, for example, 
[S]). All these make the problem for more detailed study of BH's structure in the composite modern 
non-linear theories with many interacting fields in 4D non-Euclidean space-time actual. 

A new feature of BH in these theories (even in presence of only one additional scalar field) appears 
to be the principled opportunity for appearance of several horizons of miscellaneous types with different 
space-time structure between them [5] - - [TT]. Nevertheless, there are no convincing arguments that 
one can consider the space-time domains between the various horizons as a real co-existing structures 
in a physically meaningful way, the investigation of such solutions at present becomes a part of the 
general BH problem. 

BH with three horizons arise even in the simplest generalizations of GTR - the GTR with a 
cosmological constant |12| . where historically for the first time the necessity to extend in appropriate 
way the concept of event horizon and to consider a model of space-time with complicated causal 
structure was suggested. An appearance of three horizons in space-time was obtained earlier also in 
the models of space-time with torsion 

The wide classes BH with many horizons exist in various models of the dilaton gravity (see |14j 
and numerous references there). 

BH with two event horizons were obtained numerically in the Einstein-Maxwell theory with mass- 
less dilaton ^J^j. There was pronounced conjecture, that the presence of a massive dilaton in such 
models BH with three horizons also are possible, which was numerically demonstrated in JtJ! within 
the Einstein-Born- Infold model. 

Is is well-known that "usual" BH in GTR evaporate due to quantum generation of particles by 
their strong gravitational field. This phenomenon was discovered by Hawking and provoked a rough 
progress of the quantum theory of BH. Hawking also remarked ^7], that for "exotic" BH with many 
horizons a reverse phenomenon with respect to vaporization can take place. A simple example is the 
almost singular Narai BH |18j . That is why BH with many horizons represent a special concern for 
quantum theories and were studied by many authors (see references in |14jV 

Common feature of all BH models with several horizons is their strong nonlinearity, which, as a 
rule, makes impossible their exact analytical description. As a result especially actual is the problem 
of developing an adequate numerical methods for obtaining the solutions with several horizons. Based 
on our previous studies ^UGOIj here we offer and apply such general method, which we illustrate using 
the specific BH example in the Einstein-Born- Infold model with massive dilaton. This model appears 
to be suitable polygon for development of numerical methods, which one subsequently can apply in 
other theories. In the present work, which is extension of |19|. we present the formulation and the 
numerical algorithms of solving of the problem. It is shown that the correct formulation of a BVP for 
BH equations (and consequently the performance of the method of solution) is essentially determined 
by both the number and the type of horizons. 

We solve the non-linear BVPs using iterative methods based on the continuous analog of the Newton 
method |21| in combination with method of collocation for solving of arising linearized problems. 
Compared with the methods based on solving of a Cauchy problem (see, for example, ^H]), such 
approach has definite advantages. 
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2 Formulation of the BVP 



The dimensionless BH equations have the form (the marks below are similar to those used in |19l [20 ) : 

-/' + F(r,f,<p,ip')=0, (2.1a) 
-f(<p ,, + l<p')+$(r,<p,<p') = 0. (2.1b) 

Here f(r) is the function involved in the metric space-time 

ds 2 = -f{r)e 2S{r) dt 2 + /- x (r)dr 2 + r 2 dfL, 

and <p(r) represents the dilaton field. The radial coordinate r £ [Ri, oo), where the constant Ri > 0, 
and "right hands" F and $ are set through expressions 



1 — f r 2 — , /r 4 -4- n 2 

F = - + 2 exp{2a^} - — - n 2 V(<p) - r f<p' 2 : 



ry V(v) 2exp{2aip}- 



r 



-2aexp {2o!( ( c} 



In these expressions the parameter q corresponds to the electric charge of BH, and 7 2 is a dilaton mass, 
V(<p) - potential of the dilaton field. The choice of the sign of coupling constant a = ±1 determines 
the sign of dilaton field: a = — 1 corresponds to <p(r) > 0, whereas a = 1 corresponds to ip(r) < 0. 
The case 7 = corresponds to BH with massless dilaton, investigated in works [221 EH! • Let us remark 
that in [19| the equations l|2.1(l are written in a different way by introducing an additional variable 
m(r) = r (1 — f(r)) /2 (local mass). 

For some given solution (p(r) of problem H2.1J) it is necessary to solve the following Cauchy problem 

5' + rip' 2 = 0, lim<J(r)=0, (2.2) 

r — >oo 

for the metrical function S(r). 

As usual (see, for example, pQ), positive zeroes Rh, n , n = 1,2..., Nh > 1 (Nh is a number of 
horizons) of metric function f(r) we shall call BH event horizon. We will show below that in the BH 
model under consideration, depending on values of electric charge q and mass 7, no more than Nh = 3 
horizons may exist. Greatest of them Rh we shall call an exterior (physical) horizon, and extreme left 
Ri _ Rh — an internal horizon. 

For the BVP related to equations l|2.1|l to be closed, it is necessary to formulate respective boundary 
conditions. 

First of all we shall note that for 7 ^ on the right end when r — ► 00 the asymptotic conditions 
take place: 



y/ for metric function: 
yj for dilaton field: 



2 M n 2 

/(r)->l-^ + %; (2.3) 



( 2 - 4 ) 



In formula l|2.3|l indicates the BH Keplerian-like mass. 

The posing of boundary conditions on the right end, and hence the appropriate numerical method 
for the respective BVP essentially depend on both the kind and the number of horizons. Let us consider 
in detail the basic cases of possible boundary conditions. 
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Figure 1: Metric function /(r) for several values of dilaton mass 7 



2.1 BH with regular event horizon 

We shall refer to the horizon Rh as a regular one, if the derivative 

f'(Rh)*0. (2.5) 

Let us consider the formulation of BVP for BH with single regular horizon. For Nh = 1 the semi- 
axis (0,oo) is splitted by the point r = Rh into two areas: internal Di nt = (0,Rh) and external 
D ext = (Rh, 00). In the area D ext BVP for equations 1)2.1)1 is solved using the "standard" condition 
on the horizon 

f{Rh)=0. (2.6) 

Apparently the point Rh is a point of degeneration of the Eq. 1)2. lb|) . To ensure the regularity of 
solutions on the horizon, it is necessary the equality 

$(Rh,tp h ,tf/ h ) = 0, 

to be fulfilled, which has the following detailed form 



TVfaO - ^ - 2 exp{2a^}^P±Z 



if 



2 V'(<p h )-2aexp{2a<p h }-2 v r2 h =0. (2.7) 



Here and henceforth the subscript h means that the value of the respective function is calculated at 
the point Rh- 

Let us point out that for a given mass the boundary problem Ij2.1|l - 1)2.7(1 is a problem with 
free boundary |25| . because the point Rh is a priori unknown. In many cases, however, (see discussions 
in|SJ from the computational point of view, it is more convenient to set the magnitude Rh > and 
solve the problem in D ext with a fixed boundary. Thus, mass is determined from the asymptotic 

The numerical method of solution of BVP for BH with given event horizon is presented in Section 

El 

Visual examples of metric function f(r) for Rh = 0.1 and given values of dilaton mass 7 = (curve 
noted by A), 7=1 (curve □), and 7 = 2 (curve 0) are introduced in Fig.^ Solution A corresponds 
to BH with a massless dilaton considered in articles [221 |2Sj ■ 

Let us note that in a series of works (see, for example, |22l 128) 1 to formulate the BVP for BH the 
authors confined themselves within the scope of 1)2.511 . 
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Let us consider the formulation of BVP in the internal area Di nt = (0, Rh)- Let (/+(r), ip + (r), Rh) 
be the solution of problem (|2.1|) — (|2.7|l for r > Rh and some given q and 7. Let us suppose that 
Eqs. I|2.1I) are fulfilled as well at r < Rh- To obtain the solution /_(r), ¥>-(r), i?; in some subdomain 
r G (i?;, depending on argument Ri > we require a continuity of functions f(r) and <p(r) at the 
point Rh- It leads to 

f-(Rh) = 0, (2.8) 
<P-{Rh) = f+(Rh), (2.9) 
^'_(i?,0 = <p' + {Rh). (2.10) 

For a closed boundary problem in (Ri,Rh) one more boundary condition is necessary Having in 
mind the absence of horizons for < r < Rh the choice of this condition is arbitrary enough. For 
example, such condition can be 

\f-(Rl)\ = l. (2.11) 

Since parameter Ri is a'priori unknown, problem ^2.ip . i|2.8|) - l|2.11|l appears to be a problem with 
free left boundary. Numerical algorithm for solving such BVPs is presented below in Section 0] 

2.2 BH with extremal horizons 

We shall say that BH has at point R e an extremal event horizon if the conditions 

f(R e ) = 0, (2.12) 

f'(Re)=0, 

are satisfied. 

An example of solution f(r) with extreme exterior horizon Rh w 0.776 and internal Ri = 0.1 for 
q = 1 and 7 = 3.6 is shown on Fig. ^ (curve noted by V). 

Using Eq. H2.1a[) we can write the last expression in explicit form 



1 + 2exp {2aip e }(R 2 e - t/R* + q 2 ) - R 2 el 2 V^ e ) = , (2.13) 

where subscript e means that the respective magnitude concerns the extreme horizon. 

In order to be regular on the extreme horizon the solution has also to satisfy a constrain of kind 
B2.7J1 . In this case it simply becomes 



^-y'(^ e )-2aexp{2a^} e V R2 e + g = 0. (2.14) 

For given values of physical quantities q and 7 Eqns. (|2.13|l and H2.14J) form a closed system of 
non-linear algebraic equations for determination of possible extreme horizons R e (q,j) and respective 
boundary values of dilaton tp e {Q, 7)- 

Let us note that in the presence of an additional condition for an extremeness of horizon l|2.13|l the 
system of boundary conditions is preconditioned and asymptotic condition (|2.3|l becomes "redundant" . 
From the physical point of view it means that the extreme horizon does not exist for any value of mass 

Moo- 

By exception of radical term from Eans. (|2.13|) and l|2.14|l one can obtain an explicit dependency 
between dilaton mass 7 and parameter ip e : 

Re= 1 (2.15) 

ly /V(<p e )-V>(<p e )/2a 
Obviously, potential of dilaton field V(ip) has to satisfy condition 

v(< Pe )-^ e v'(i Pe )>o. 
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Figure 2: Functions C(tp),C\{tp) and Ci(tp) for 7 = 1. 



In particular, with a = — 1, and 

V(tp) = tp 2 , (2-16) 

the last relation becomes 

+ <Pe > 0. 

In this manner horizon R e is determined for tp e € (0,oo) and 

Re= r^r— - (2-17) 
To compute the unknown parameter ip e one should use Eq. 



7 2 1 
1 - -j- cx P{-2a<^ e } — V'(tp e 
4 2a 



+ 



q 7 exp{2a<,9 e } 



la 



0, (2.18) 



following from dependencies (|2.14(l and l|2.15|l . 

Simple notion about the qualitative behavior of solutions l|2.18Jl can be obtained in particular case of 
quadratic dilaton potential l|2.16(l . Then Eq. I|2.18|l can be written in the form (for simplicity subscript 
e is omitted): 

C(tp, q, 7) = Ciiy, 7) - C 2 (tp, q, 7) = 0, (2.19) 



where 



Ci(p,7) = 1 + —ipexp{2(p}, C 2 ((p,q,j) = exp{-2</?} tp (1 + if) 2 . 



For given electric charge q and dilaton mass 7 the function Ci(tp) monotonically increases when 

tp E [0,oo), and C(0) = 1 and lim Ci{tp) = 00. By contrast, the function C^fy) monotonically 

tp — ► 00 

increases in interval [0,1) and decreases in (l,oo). Since 6*2(0) = 0, C 2 (l) = and C^l) < 0, for 

tp = 1 the function C^fy) has a maximum. Because of lim C2[tp) = 0, Eq. (|2.19|) can have on interval 

— >oo 

(0,oo) no more than two roots, the greater of which we shall denote tp £i i, and the smaller — <p>e.i 
(hereinafter, if it does not result in misunderstanding, we shall suppose tp e ,\ = tpi and tp e ^ = tfh)- 
These reasonings are illustrated graphically in Fig. where the dashed line indicates the plot of 
function C\{tp), the dash-and-dot lines indicate the plots of function C2(tp) for three values q = 1.6, 
q = 1.79 and q = 2.3, and the continuous lines — the graphs of respective dependencies C\(tp). In 
accordance with formula H2.17fl the smaller root tp h < 1 of Eq. (|2.19l) corresponds to exterior extreme 
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Figure 3: Relationship R (q) 



Figure 4: Relationship R (7) 



horizon Rh of BH and the greater root ip e — to the BH with internal extreme horizon (see below the 
formulation of respective BVP). In order to magnitude p e > 1 (the root is located to the right of the 
maximum of function C2 ((/?)) the following condition should be satisfied 



where e ~ 2.718 . . . In particular, for 7 = 1 we obtain q > 2.293 (see Fig.EJ. 

For given values q and 7 it is easy to solve numerically the non-linear algebraic equation l|2.19|l 
and compute the values of the roots <fu{<lil) and ^(9,7). As an example on Fig. [21 the dependencies 
Rh (q) (continuous curves) and Ri (q) are shown (dotted curves) for two values of dilaton mass 7 = 1 
and 7 = 2. Obviously when charge q of exterior extreme horizon increases, Rh increases linearly with 
coefficient ~ 1/2, and internal Ri decreases. The points indicated through D correspond to the triply 
degenerated horizons Rh = Ri = Rd- 

Similarly on Fig. 0] the dependencies Rh("f) (continuous curves) and Ri{^) (dotted curves) for two 
values of a charge q = 1 and q = 2 are demonstrated. It is seen that the main variation both exterior, 
and internal horizons under influence of dilaton mass 7 is localized in some neighborhood of the point 
of triply degeneration D. The large values of dilaton mass 7 render minor influence on the magnitude 
of horizons. 

2.2.1. BH with exterior extreme horizon 

Let us suppose, that for some q and 7 Eq. I|2.19|l has two roots and consider the smaller of them 
iphiq,^), to which according to formula l|2.17[) corresponds an external horizon Rh- The setting of 
magnitude Rh actually means that the left boundary of the area D ext is known and, therefore, system 
(|2.1|) on interval D ext is solved with boundary conditions (|2.6(l and 



on the horizon, and also with dilaton asymptotic (|2.4J) on the right end. The BH mass is obtained 
from the asymptotic of metric function l|2.3|l . 

An example of a solution obtained numerically with external extreme horizon is shown in Fig. ^ 
(see the curve noted by V). In point E the plot of the function /(r) concerns a horizontal tangent, 
i.e., the condition f'(R e ) = is fulfilled. 

For limited values of dilaton mass 7 BH has also a regular internal horizon (see below the discussion 
of results of the numerical experiment in Section [SJ). This horizon Ri is an unknown, and, hence, in 
area D m id = (Ri,Rh) it is necessary to consider a problem with free left end for Eq. I|2.1|l . Let 
{/(r), ip(r), Rh} be solution of the problem in the exterior area D ext . in addition we suppose function 
/(r) to be continuous, and function <p(r) is smooth at point Rh- Then for Eq. H2.1(l on interval D m id 




<p(R h ) = fh(q,l) 



(2.20) 
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Figure 6: ip — discriminant. 



the boundary conditions look like 1)2.6(1 . ((2.20(1 at Rh, and also ((2.6(1 and 1(2.7(1 at Ri (in the last two 
expressions one should substitute for subscript h with 

2.2.2. BH with internal extreme horizon 

As was mentioned above in Section ^. 2l the BH extreme horizon Ri can be internal, e.g., R &i i = Ri < Rh- 
Let us consider the formulation of BVP for Eq. 1(2.1(1 in this case. 

Points Ri and Rh divide the half curve into three areas: internal Di nt = (0, R{), intermediate 
D m id = (Ri,Rh) and exterior D ext = (Rh, oo). Let us at first consider the problem in interval D m id- 
At point R e along with the two conditions of kind ((2.121) we impose an additional one 

l p(R l ) = < P l(q n ), (2.21) 

where the value <p e (q, l) is the greater root of Eq. 1(2.19(1 . At the unknown right boundary Rh both the 
condition of existence of horizon 1(2. 6|) and regularity condition 1(2.7(1 should be held. Thus the problem 
in the area D m id is self-contained. 

Let the solution {/(r), ip{r), Rh} in this area is found. Then, assuming the solutions at point Rh 
to be continuous in area D ext , it is necessary to solve Eqns. (12.1(1 with boundary conditions (12.6(1 and 
(|2.20(l on the left end, and also with 1(2.4(1 on the right one. 

Particular example of numerical solution f(r) with internal extreme horizon Ri as 0.84 at the point 
E for charge q — 2 and dilaton mass 7 = 1 is presented in Fig. |S| 

The formulation for seeking the solution in the internal area Di nt is similar to depicted above at 
the end of Subsection 12. II 

2.2.3. BH with triply degenerated horizon 

Let for some values of charge q and dilaton mass 7 Eq. ((2.19|) has two roots iph < <pi (see Fig. EJ) - As for 
ip G (iph,ipi) the inequality C\((p) < C2 (^) is satisfied, in some point ip m (and in view of a continuity 
and in some area containing point ip rn ) in this interval function C(ip) has a negative minimum. When 
dilaton mass 7 decreases, the graph of function C(<^, 7) rises up, the graph of function Ci (1/5,7) falls 
down, and the distance between roots <ph and <p>i decreases. For some critical 7 = jd the minimal point 
(fid = fm of function C(y>, 7) concerns to horizontal axis being a single root of Eq. 1(2.1911 . Thus the 
relations 



C{<fd,qd,ld) = 0, 
dC , 



dip 



(2.22a) 
(2.22b) 



take place. 
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Figure 7: An example of a solution f(r) with triply degenerated horizon. 



Eqns. dete rminc on plain (5,7) a smooth curve (see Fig. EJ|, which is called (^-discriminant 

of the function C(<p,q,j). The equations of discriminant in parametric form can be obtained from 
definition lf!T52j) 

Id = CXP( ^V 2(1-^), (2.23a) 

'fid 

_ cxp(2y d ) 

qd 0/1 2 - (2.23b) 



The simple calculations demonstrate that relation 12.22bjl expresses the condition of a vanishing of 
second derivative of metric function f(r) in horizon Rd 

f"(R d ) = 0. 

Such horizon is called in articles ^3 EH| triply degenerated. 

Thus, the algorithm for solving the problem for BH with triply degenerated horizon looks like the 
following. It is more convenient to set dilaton mass 7 = 7<j. Then according to formula (|2.23a|l it 
is possible to compute < (fdd) < 1, and by means of relationship l|2.17|l - the horizon Rd of BH. 
The BH charge q is found from i|2.23bjl . Therefore, in external D ext /internal Di nt domain BVP for 
triply degenerated horizon seems to be a problem with fixed left/right boundary Rd, where one sets 
conditions of kind Q2.6|l and H2.20JI . If the solution of this problem is found, then respective BH mass 
Moo can be found from the asymptotic expression (|2.3(l . 

An example of triply degenerated solution f(r) obtained numerically for qd ~ 0.836 and 7^ = 4 
is presented in Fig. Point D (Rd ~ 0.456) is an inflection point for metric function /(r), i.e., in 
addition with f(Rd) = the conditions f'(Rd) — and f"(Rd) — are fulfilled also. 



3 Numerical Method for Solving a Problem with Fixed Bound- 
ary 

As was mentioned above, in a case of BH with exterior extreme horizon Rh it is necessary to decide a 
BVP for Eq. Ij2.1|l in domain D ext with fixed physical parameters at the right end Rh and system with 
three boundary conditions H2.6(l . (|2.20(l and (|2.4|l . For the numerical implementation it is expedient 
to apply the continuous analogue of Newton method (CANM) which is very suitable for treating of 
numerous problems in physics |21|. 

Let us introduce an "actual" infinity i?oo < 00 and consider the non-linear functional equation 

x(y) = o (3-i) 
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with respect to pair y = {f,<p}, y £ C 1 [Rh, Roo] x C 2 [Rh, Roo], where the coordinates X™(y), n = 
1, 2, . . . , 5, are determined as follows: 

X (1) = -f' + F(r,f,<p,<p'), (3.2a) 

X (2) = -/fv"- W *(»•,¥», V'), (3.2b) 



X (3) = Wh), (3.2c) 
X (4) = <p(R h ) - Vhfa'r), (3-2d) 
X^ = + (3.2e) 

Let Eq. 1[1 have an insulated solution y* = {/*(r), </?* (r)} and initial approximation j/° = 
{/ (r),y (f)} near to solution y* is known. For small values of charge q and dilaton mass 7 we 
can use as initial approximation asymptotics 12.. '{[I and i|2.4|) . 

Let us introduce continuous parameter t £ [0, 00) and assume trajectory y(t) satisfying the abstract 
Cauchy problem 

X'(y(t))^-+X(y(t)) = 0, 2/(0) = 2/0 . (3.3) 
at 

Here x'{v) is Frechet's derivative of the non-linear operator x{u)- -"- n |24j it is shown that when the 
operator x{y) IS smooth in some vicinity of the sought solution y* the limit relationship 

lim - 2/*|| ►O. 

t — >oo 

takes place. 

For the numerical solution of Cauchy problem l|3.3|l it is easy to take advantage of the explicit Euler 

method on non-uniform grid t k +i — t k + r k with variable step Tfc, k = 0, 1, . . . As a result we come to 
iterative process 

x'{yk)w k = -x(Vk), (3.4) 

Vk+i = yk+T k w k , (3.5) 

which allows on each iteration k by help (|3.4J) to compute next approximation y k +i to the exact 
solution. When r k = 1 we come to the classic Newton method. 

In the convergence of iterations (|3.4|) . (|3.5|l to the solution of the evolution Cauchy problem 
(|3.3II in finite time interval t when r k — > is proved, if in vicinity of exact solution y* operator x(j/) 
satisfies some additional conditions. 

The variation of step T k can be ruled during the iteration process |21l I2H| . 

On each iteration Eq. 1)3.4(1 in the problem under consideration is equivalent to the following linear 
problem for the coordinates of correction vector w(r) = f](r)} (for simplicity hereinafter index k 

will be omitted): 

dF dF dF 

^' + a^ v ' + W^ + *p> f " F (r ' f,lf,tp ] ' (3 ' 6a) 
( „ i A 0* , / „ 1 ,v 0$ 



, , (3-6b) 
f[<p" + -<p')-*(r,<P,<p') 

£(Rh) = -f(Rh), (3.6c) 

J]{Rh) = <Ph(q, 7) ~ <p{Rh), (3-6d) 

2 

ryORoc) = -<p(Roo) ~ -gr- (3.6e) 
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It is convenient to select initial approximation {fo{r),ipo(r)} satisfying initial conditions 1)2. 4[1 . 12.6)1 
and 1)2.20)1 . Then on each iteration right hand sides l)3.fic)l - ()3.fie)l are equal to zero, and the contribution 
of boundary conditions in error <5(t^) will be zero, respectively. 

Let us suppose that {C( r )i 7 7( r )} is a solution of the problem (|3.6(l . Then the next approximation 
to the exact solution is computed from relation (see formula (^3.5|) ^ 

f k+i = f k +r k e, v k+i = <p k +T kV k . (3.7) 

In cases, when horizon Rh < R e ,i or Rh > R e ,2, dependence i?/ l (M co ) appears to be univalent (see 
below the discussion in Section EJ). This allows to consider the quantity Rh as a parameter and thus 
substitute the problem with free left boundary problem with fixed boundary Rh. In this case the right 
hand side of the expression (|3.2d|l becomes form similar to Eq. 1)2.7)1 : 

X W = *(Rh,<Ph,<p' h ) 



and boundary conditions 1)3. 6d|) is substituted by 

r)'{R h ) -$(R h ,9h,Vh)- 



<9$ 

dip 



<9$ 

r]{Rh) + 5-7 



I!,, 



On each iteration k a linear BVP (|3.6)l is solved in finite interval (Rh, -Roo)- For discretization the 
method of collocation at the Gaussian knots of grid exponentially condensed to horizon Rh is used. 

Let u* be the exact solution of the original continuous problem l|3.1|l on a finite interval of time, 
u* h — the exact solution corresponding to the discretized non-linear problem 

Xh{u h ) = (3.8) 

in a finite interval of time, u k - approximation to Uh after fcth iteration when condition < £, 

where < e <C 1, is fulfilled. Let us note that if the discretization method does not vary from iteration 
to iteration, then the grid representations for 1)3.6)1 follow from grid representation (|3.8|) of original 
equation l|3.1|) . 

For convergence estimate of the method we shall consider the inequality 

||y* - ul\\ h < \\y* - u*\\ h + \\u* - u h \\h + IK - u k h \\ h ■ 

It is possible to show j25 that ||m* - u h \\h < 0(h r ), \\u* h - u^W < Bh\\Xh(Uh)\\> w ^ ere B h - some 
constant. Then for Bh \Xh{u^)\\ <C 0(h r ), which is fulfilled for enough small e, the precision of obtained 
approximate solution is specified by two first summands on the right hand side in the inequality. 

The error 5^ — \\y* — u*\\h is investigated numerically on a fixed grid for different values of actual 
infinity i?oo, which one select, so that <5oo remained small with respect to two other summands. Thus, 
the precision of the approximated solution is close to the theoretical estimate of method of difference 
approximation of Eq. 1)3.1(1 . 



4 Numerical Method for Free-Boundary Problem 

An exhausting enough review of methods for solving free-boundary problems is presented in treatise 
|25| . In present work the boundary- value problem for Eqns. 1)2.1)1 with free external Rh and given 
extremal internal Ri horizons is rendered to non-linear eigenvalue problem with spectral parameter 
Rh, which in turn is solved by CANM. Let us remark that such approach has good reputation for 
treating various problems in astrophysics [261 127| . the theory of Josephson junctions 28. , etc. 

A formal shortcoming for solving the problem with unknown external horizon by CANM is the 
absence of explicit dependence between the equations and boundary conditions and horizon Rh- In 
order to enter explicitly parameter Rh we introduce new variable x under the formula 

Uh — rLi 
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After the change of variables (|4.1|l interval [Ri, R^] renders to [0, 1], at that d/dr = (Rh~Ri) 1 d/dx. 
Let us assume z = {y,R h }, y = {f,cp}, z € C^R^R^] x C 2 [R h , R^} x K. Then BVP for Eq. flTT) 
can be written similarly to Eq. I|3.1[l : 

x(y,Rh) = 0, (4.2a) 
N(y) = 0, (4.2b) 

where vector x( z ) is specified using expressions (to avoid the introducing of redundand notations 
further in this section we substitute (.)' = d(.)/dx) 

X (1) = -f' + F[xJ,<p,tp',R h ]=0, (4.3a) 

x (2) = -/ (V + -<p')+*[x,<p,<p',R h )] = o, (4.3b) 



x 

.(3) 



X w ^ /(0), (4.3c) 
X (4) = /(I), (4.3d) 
X ^ = ${l, Vh , V ' h ,R h ), (4.3e) 



and the left part of "norm condition" l|4.2b|l has the form 

N(y) = v(0)-Vi(q,j)- (4.4) 

Here through F and i> the right hands of the BH Eqs. 12.1f> after substitution in them l|4.1|l are 
indicated: 

F = {Rh-R^F^ix.R^J^^Rh-R^- 1 ^ 1 }, 
$ = {Rh-Rif^^ix.R^^^Rh-R^- 1 ^'} . 

As in Section|31 we introduce continuous parameter t S [0, oo), supposing the CANM relation being 



x' y (y,Rh)w + x' Rh (y,Rh)p + x(y,Rh) = 0, (4.5a) 

N'(y)w + N(y) = 0, (4.5b) 
w = y, p = R h . (4.5c) 

Problem (|4.5a|) - 1)4. 5c|) has to solved with initial conditions 

y(0) = y , R h (0) = R h , O - (4.6) 

The numerical realization of CANM can be implemented by different methods for approximated 
solving of differential equations. Let us consider further Euler iteration process corresponding to 
Cauchy problem l|4.5f) on non- regular generally grid ifc+i = tfc + Tfc, fc = 0, 1, . . At fc-th iterations 
it is necessary to solve the linear operator equations (|4.5a() and 14.5bl) . after that the subsequent 
approximation to the exact solution is found through the formulas, following from relations (|4.5cll : 

= /' + Tk w k , R k h +l = R k h + r k p k . (4.7) 

We shall search the solution of linear equation l|4.5a|l in the form (for simplicity henceforth we omit 
iteration index k) 

w = u + pv , 

where u(x) and v(x) are new unknown functions. Substituting this decomposition into 1)4. 5a|) and 
equating coefficients we receive the system 

x' y (y,Rh)u = -x(y,Rh), (4.8a) 
x ' y (y,Rh)v = -x' Rh (y,Rh). (4.8b) 
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Figure 8: Dependence Rh(Moo) for small 7. 



Figure 9: Dependence i?/ i (M tx) ). 



Let the solutions of these equations be computed. Then the derivative p can be obtained from 

p = -[N'(y)v}- 1 [N(y) + N'(y)u}. 



(4.9) 



It is expedient to select the initial approximation yo so that it obeys norm condition l|4.5b|) . In that 
case expression (|4.9I) is simplified: 



P =-[N'(y)v] X N'(y)u. 



(4.10) 



Let us give the explicit form of Eas. l|4.8|) and (|4.9|l . Let u = {uj(x), u v (x)}, and v = {vf(x) 7 v ip (x)}. 
Then we have 



dF dF dF 

-u' f + - — ' + — tl/ + — — u v = f — F(x, /, ip, ip , Rh) , 



; dp' v ' df~ J ' dp' 
„ 1 



<9$ 

dp 



\ <9<i 




1 + w 




f 


(, 










+ a^ 7 


x= 



u'Jl)=-$(l,<p,<p' t R h ), 



(4.11a) 

(4.11b) 

(4.11c) 
(4.11d) 

(4.11e) 



dF 



dF dF dF 
~df Vf + 9V VV ~ ORii 



ft " jl 1 A , a * ' f " , 1 A , 5 * 



dp 



«/(0) = 
v f (l) = 
d$ 



v'M) 



9$ 



9Rh 

yfar) - y(Q) - ~M°) 

M0) 



(4.12a) 

(4.12b) 

(4.12c) 
(4.12d) 

(4.12e) 
(4.13) 
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Thus, the iteration process is conducted in following order. For given initial approximation fo(x), 
tfo{x) and Rh we compute increments u(x), v(x) solving boundary problems (|4.11|) and l|4.12[l . Further 
by formula l|4.13[) we hnd increment p of horizon Rh . The next approximation of the exact solution is 
obtained using (|4.7J) . 

The linear BVPs (|4. llf> and (|4.12(l are solved numerically through the collocation method of order 
0(h A ) in the Gaussian knots of grid, which is exponentially condensed to horizon Rh [2E]- Let us note 
that the left-hand sides both the equations and the boundary conditions are identical, what simplifies 
the solution of respective matrix problem. 

Let the solution in domain D m id be found. Then supposing a continuity of functions at point Rh, 
the problem in external domain D ext is treated as a problem with fixed left boundary (SectionEJ. If it 
is necessary to find the solution in internal domain Di nt (see the end of Section^ the corresponding 
BVP with the free left boundary can be solved as is depicted above. 

Apparently the solution of the problem with extremal external and regular internal horizons can 
be obtained by methods presented in Sections and 0] 

5 Discussion of Numerical Results 

^From physical point of view very relevant is the link between BH mass and horizon Rh for 
different values of parameters q and 7. 

The numerically obtained dependence Rh{M 00) for small dilaton mass 7 = 0.01 and q < 1 is 
demonstrated in Fig. El F° r such values of 7 it is necessary to expect that uniqueness of the dependence 
R(M OQ ) and linearity for large values are saved in a wide range of variation of charge q. This 
deduction is similar to the recent results |22l |2~3*] concerning the case of massless dilaton. It is easy 
to calculate that the triple degeneration of horizon, which initiates an essential influence of mass 7 on 
number and form of horizons (see below), takes place for q > 130. In this sense the curves in Fig. OH 
have test nature. 

The influence of finite dilaton mass 7 is shown in Fig.Elfor given charge q = 1. When 7 < 7^ ss 2.65 
(see curves indicated through and □) dependencies R(M oa ), though hardly distorted, remain single- 
valued. Curve R{M oa ) (denoted through O) corresponding to extremal value jd, has in the point of 
triple degeneration of horizon D (Rd « 0.58) vertical tangent (for inverse function M^Rh) point D is 
an inflection point). The corresponding solution for metric function /(r) is plotted in Fig. Q curve □. 

For further increase of dilaton mass 7 > jd BH has (see Section l2~2")l two extremal horizons Ri(q, 7) 
and Rh(q,j), Rh > Ri- In the graph of dependence R(M OQ ) the presence of horizons Ri and Rh 
expresses in appearance of typical ^-shaped curves (in given case indicated through V and A). Two 
different solutions of Eq. (|2.1() with different masses M^q,^) and M/(g,7) correspond to these horizons. 
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Figure 11: An example of solution f(r) with Figure 12: Examples of solutions f(r) 

two regular horizons. with three regular horizons. 



Let us fix the value 7 = 3.2 > 7^ (curve noted by symbol V in Fig.EJ). In Fig.E3this curve together 
with curve A are plotted in an appropriate scale. Moving on arc ABE1E2C from left to right the 
dynamics of variation of number of horizons apparently looks like the following. From point A up to 
point B± BH has single regular horizon Rh, on vertical straight BE2 - two horizons, and horizon in 
point B is regular, and in point E2 - external extremal one ( Rh ~ 0.72). Further, for given M^ on 
the segment between vertical straights BE2 and E\C BH has three regular horizons, on the straight 
line E\C - two, and vertical coordinate Ri « 0.42 of point E\ is an internal extremal horizon. At last, 
on the right of point C BH has again only one regular horizon (R = Rh). 

The "motion" of extremal horizons Ri and Rh as a result of the variation of mass 7 for values of 
electric charge q = 1 and q — 2 is presented in Fig.^J Similar dependence of the magnitude of horizons 
from charge q for fixed 7 = 1 and 7 = 2 is demonstrated in Fig. 01 

For large enough 7 the vertical coordinate Rh of point E2 increases linearly (see Fig. 0J, and the 
coordinate Ri of point E\ comes downward. Thus the point Mo (mass specified by the electrical field 
and the dilaton), fulfilling the formal equality R(Mo) ~ 0, can be located on the right to the vertical 
straight line BE2, i.e., will be executed Mo > Mh (see the curve noted by A in Fig. 110(1. This means, 
that for respective 7 (in particular case 7 = 5) the distribution of horizons is as follows: at point E2 
BH has a single extremal horizon (see Section l2~2"|) . Further, in the domain between vertical straights 
through E2 and Mo BH has two regular horizons. In the right of the point Ma up to E\ BH has three 
regular horizons, at point E\ - two horizons Ri and R c = Rh, and internal Ri is extremal. At last, for 
large enough masses M x (in the right of point E\) BH has a single regular horizon. 

In Fig. II li the examples of solutions with two regular horizons in neighborhood of extremal solution 
A, Moo ~ 3.96 (curves noted by characters □, and V) are plotted. At that computation values of 
the BH mass: □ - Moo = 4.1, - Moo = 4.22 and V - M^ = 4.36. 

In Fig. 1121 two examples of solutions with three horizons for 7 = 1, q = 1.9, Rh ~ 1-48 (the dotted 
curve noted by a character V) and q — 2, Rh ~ 1.69 (continuous curve A) are shown. 

It is easy to obtain a condition, for which on the leftward of vertical straight line BE2 through point 
Rh (see Fig. BH does not have any horizons. For this purpose we mark by M; and Mh the masses 
of two BH, for which Ri and Rh are respectively internal and external extremal horizons. Let us build 
Hcrmitian cubic polynomial S3 (R) interpolated the function M^ (R) on the interval R S [Ri, Rh] ■ 
Such a polynomial is based on known conditions M oc (i?;) = M;, M ao (Rh) = Mh, M '^(Ri) = 0, 
M '^(Rh) = 0. Let us also set t= (R- R t )/AR, where AR = R h - Ri and AM = M\- M h . Then 
the polynomial S3 (R) becomes 

S 3 (R) = AM (2t 3 - it 2 ) + Mi, iG[0,l]. 

For further purposes it is convenient to rewrite the equation for local coordinates t corresponding 
to some given mass in the form 

2t 3 -M 2 +fi = Q, (5.1) 
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Figure 14: BH temperature. 



where the reduced mass p is evaluated under the formula p = (Mi — M DO )/AM. The form and values 
of the roots of this equation depend on the magnitude of p, i.e., from the relation between three masses 
Moo, Mi and Mh- It is not difficult to verify that when < p < 1 the equation (|5.1|l has three real 
roots placed on the interval [—1/2, 3/2]. 

As it is seen from Fig.^3 for large enough 7 BH can have a regular horizon, only if mass Mq < Mh, 
which is equivalent to ji < 1. If ji > 1, then according to the above-stated remarks BH has two or 
three horizons. We continue formally polynomial Ss(R) on the segment R £ [0, R{\ and require mass 
Mq corresponding to local coordinate to — —Ri/AR, to coincide with extremal mass Mh- In this case 
coefficient p = 1 and Eq. (|5.1|l possesses one simple root to = —1/2, and also two repeated roots t = 1. 
Further we introduce into consideration quotient p = Rh/Ri > 1- Value p = 1 corresponds to the case 
of BH with triply degenerated horizon, i.e., to the point D in Fig. EI For 7 > 7<j the inequality p > 1 
holds. In particular, for to the magnitude p possesses the value 



Rh{q,l) 



= 3. 



The respective dependence between charge q and dilaton mass 7 is illustrated in Fig. 

In this way if 1 < p < 3, then to the left of extremal mass Mh there is a single- valued branch of 
dependence Rh(M OQ ). In interval Mh < Mo < Mi (point Mo is located to the right of point Mh and 
p > 3) parameter p £ (0, 1). From here it follows t £ (—1/2, 0) and, hence, BH has even two regular 
horizons. 

For calculus of coordinate R c of point C we continue polynomial S3 (R) on the segment R € [Rh , R c ] > 
i.e., t e [1, (R c - Ri)/AR}. Then M c = Mi, p = and Eq. (|5T|) has single nonzero root t c = 3/2, to 
which corresponds horizon 

Rc - ■ (5 - 2) 

The appearance of a triply degenerated horizon with consequent arising of extremal horizons under 
influence of dilaton mass 7 one can see in the graph of dependence between BH temperature 1 



T{R h ) 



47T 



exp{-5 h }f' h 



and external horizon Rh when charge q is fixed (See Fig. 114(1 . If BH mass 7 < 7<j, the graphs of 
temperature T(Rh) are smooth curves. For 7 = 7^ the graph T(Rh) concerns the abscissas in point 
Rd, and the BH temperature is equal to zero. If BH mass 7 > 7^, then the temperature curve is 
splitted to two branches corresponding to arcs AB (left branch) and E^C (right Branch) in Fig. 1101 

1 Function <5(r) is a solution of the Cauchy problem 12.21 
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6 Conclusion 

In this paper we show, that depending on the form and the number of horizons for equations of 
symmetrical charged BH in the frameworks of the string Einstein-Born-Infeld model with massive 
dilaton it is necessary to state different multipoint BVPs with unknown boundaries. For their solution 
we propose an effective iteration method based on CANM in combination to the method of collocation 
for treating the arising linearized problems. 

The horizons for the extremal BH solutions are found from non-linear algebraic system depending 
on charge q and dilaton mass 7. On the graph of horizons against BH mass the presence of extremal 
horizons is exhibited in appearance of separate branches corresponding to BH with one, two and 
three horizons. The relevant special case is the single solution of this system, representing the triply 
degenerated horizon. After solving BVP for BH with extremal horizons and calculation of respective 
masses of BH the Hermite cubic polynomial is built. The real roots of this polynomial specify the 
number and form of horizons in general case. 

The relation of the above described BH solutions with the physical reality remains an open problem. 

After the Russian version of this talk was published in |2U], an e-print by T. Tamaki was published 
- |30| . It overlaps some part of our work. 
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